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ABSTRACT 

Challenges to computational aerothermodynamic ( CA) simulation and validation of hypersonic flow over plan- 
etary entry vehicles are discussed. Entry, descent, and landing (EDL) of high mass to Mars is a significant 
driver of new simulation requirements. These requirements include simulation of large deployable, flexible 
structures and interactions with reaction control system (RCS) and retro-thruster jets. Simulation of radiation 
and ablation coupled to the flow solver continues to be a high priority for planetary entry analyses, especially 
for return to Earth and outer planet missions. Three research areas addressing these challenges are empha- 
sized. The first addresses the need to obtain accurate heating on unstructured tetrahedral grid systems to take 
advantage of flexibility in grid generation and grid adaptation. A multi-dimensional inviscidflux reconstruction 
algorithm is defined that is oriented with local flow topology as opposed to grid. The second addresses coupling 
of radiation and ablation to the hypersonic flow solver — flight- and ground-based data are used to provide 
limited validation of these multi-physcis simulations. The third addresses the challenges of retro-propulsion 
simulation and the criticality of grid adaptation in this application. The evolution of C A to become a tool for 
innovation of EDL systems requires a successful resolution of these challenges. 

1.0 INTRODUCTION 

The ability to provide simulations of hypersonic flows over planetary entry vehicles with well-defined uncer- 
tainties on the predicted aerothermal environments is critical to the success of solar system exploration. As 
mission needs evolve the demands on simulation of hypersonic flow over planetary entry vehicles become more 
complex and the validation of those simulations becomes more challenging. Consider the role of computational 
aerothermodynamic s (CA) in the definition of aerothermal environments for representative entry vehicles in 
recent years[l-3] and the role it must assume in the near future. 

The Mars Pathfinder program (entry 1996) was the first to use numerical simulation as the primary source 
for defining the aerothermal environment for entry to Mars [4]. The vehicle forebody geometry was a scaled 
Viking shape. Viking entered the Martian atmosphere from orbit with an offset center of gravity and used a 
reaction control system to trim at 1 1 degrees angle of attack. In contrast, Pathfinder used a spin stabilized, 
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direct entry at zero degrees angle of attack. The Pathfinder simulations were validated with ballistic range 
data and wind-tunnel data obtained for Viking. The pre-flight discovery of a static instability associated with 
sonic line movement (a condition never identified in the Viking experimental database) provided a serendipitous 
opportunity to validate the prediction in flight through accelerometer data. The aeroshell was not instrumented 
but the accelerometers recorded vehicle oscillations as the sonic line jumped from the corner to the nose and 
then later again as the sonic line jumped back from the nose to the shoulder. These jumps, associated with 
evolving gas chemistry, were clearly defined within the matrix of simulations performed to define the aerodata 
book. 

Support of aerothermal design for the Orion entry capsule and more recent robotic entry vehicles to Mars 
have brought new simulation and validation needs to the forefront. Reaction control system jets fired from 
the base of the vehicles for active control authority interact with the high speed flow expanding around the 
corner of the vehicle into the wake. [5] Simulations are executed to quantify the magnitude of jet interaction 
effects on aerodynamic coefficients and to quantify the magnitude of any hot spots on the aft heat shield that 
may result from these interactions. At present, such simulations are not validated. These jet interactions are 
typically unsteady. The ability to adapt the grid to the plume shocks and free shear layers is limited by the use of 
structured grid topologies. The simulations are valuable in that model perturbations can be executed to provide 
a better understanding of the sensitivities introduced by the interactions. But clearly, more advances are needed 
if CA simulations are to independently address these challenges with well defined prediction uncertainties. 
Another issue in current design efforts is the prediction of transitional and turbulent interactions associated 
with tension ties and compression pads in the Orion heat shield. [6, 7] These perturbations to a smooth thermal 
protection system (TPS) introduce a transitional or turbulent wedge of flow downstream. Transition criteria 
are defined from ground-based tests[8] with some corroboration from a shuttle flight experiment. [9] Algebraic 
turbulence models still do the best at predicting heating levels for attached, fully developed turbulent flow but 
these methods are not suited for simulating the transitional / turbulent wedge behind such perturbations. 

Aerothermal design support for future vehicles must address new challenges introduced by: (1) entry, de- 
scent and landing (EDL) for heavy mass to Mars (one metric ton is the approximate limit for existing archi- 
tectures); and (2) limitations on payload mass to the outer planets or return from Mars which are derived from 
extreme heating levels. As will be discussed in the next section, both of these challenges are related to a desire 
to minimize system mass by utilizing aerodynamic drag to slow the vehicle. In both cases, packaging more 
mass into vehicle diameters constrained by present launch systems tends to increase the ballistic coefficient of 
the entry vehicle. In the case of large mass to Mars, if the ballistic coefficient is too large the entry vehicle runs 
out of altitude before it runs out of velocity. In the case of high entry velocity (to outer planets or return to 
Earth), a high ballistic coefficient tends to increase peak heating rates and requisite thermal protection system 
mass. The ballistic coefficient, (3 = mj (C^A) where m is mass, Cd is drag coefficient and A is reference area, 
can be decreased by increasing the reference area with a deployable structure (semi-rigid or inflatable). Several 
scenarios for deployment in the hypersonic and supersonic domain have been suggested. Simulation capability 
must evolve to handle both small and large scale (possibly unsteady) structural deformations in defining the 
aerothermal environment. Another option is to carry some retro-propulsion capability but search for optimized 
configurations that minimize propulsion system mass (and possibly thermal protection system mass) through 
engagement in the supersonic (and possibly hypersonic) domains. 

In the view forward, it is clear that CA must be able to simulate a variety of flow physics phenomena coupled 
with structural and material response to the aerothermal environment. Furthermore, the simulation setup should 
be sufficiently simple and robust to accommodate adjoint-driven optimization studies so that EDL architecture 
design trades may be implemented using best-available, high-fidelity analyses. The simulations must resolve 
jet interactions with flow expanding into the wake (RCS interactions) and opposing flow on the windward side 
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(retro-propulsion). They must also accommodate large scale surface deformation (inflatable structures). These 
topologically complex interactions will require robust grid adaptation. The simulations must also accommodate 
unsteady flows associated with these interactions. We frame these challenges in the context of simulation for 
innovation. Simulation for innovation has three defining features: (1) conceptual systems can be quickly set 
up for simulation and optimization, (2) robust grid adaptation to a user specified discretization tolerance is 
available, and (3) credible estimates of simulated aerothermal environment uncertainties are published. Feature 
(1) is a requirement shared across all speed domains but features (2) and (3) have challenges unique to the 
hypersonic flight domain. In the material that follows we will review recent efforts to address the numerical 
and physical modeling challenges required for CA simulation for innovation of planetary entry vehicles and we 
will touch base with existing flight and ground-based data to validate these simulations. 

2.0 PROBLEM DEFINITION 

The suite of physical models required to define the aerothermal environment of an entry vehicle are a function 
of relative velocity at entry interface, vehicle size and ballistic coefficient 0 = mj (CdA) (A cannonball has 
a very high ballistic coefficient, a beach ball has a very low ballistic coefficient.) In general, planetary entry 
vehicle simulations must include the capability to model gas- surface interactions (catalysis of atomic species 
to molecular species), ablation (thermal protection system (TPS) mass emanating from the boundary due to 
material degradation at high temperatures), and radiation (energy flux originating from the high temperature 
shock layer). If the vehicle is sufficiently large, the aerothermal environment at peak heating may be augmented 
by a turbulent boundary layer. In addition to these processes long associated with the simulation of aerothermal 
environments, the demands of entry, descent, and landing at Mars have introduced new challenges. (At least 
new to EDL - the challenges discussed below are also relevant to innovative concepts for aerocapture.[10]) In 
this chapter we review the driving issues that define the additional suite of simulation capabilities required for 
computational aerothermodynamic (CA) analyses of planetary entry vehicles. 

2.1 Entry, Descent, and Landing — EDL 

Braun and Manning[l 1] present EDL challenges at Mars and the need to expand beyond heritage systems. They 
note that there have been five successful landings on Mars including Vikings I and II (1976), Mars Pathfinder 
(1997), and Mars Exploration Rovers Spirit and Opportunity (2004). All five of these successful systems 
have landed masses of less than 0.6 metric tons. They have landed at low elevation sites (below 1 km Mars 
Orbiter Laser Altimeter (MOLA)). All accepted a relatively large uncertainty in landing location of 100s of 
kilometers. All of the current Mars missions have relied on large technology investments made in the late 
1960s and early 1970s as part of the Viking Program. They utilize the aerodynamic characterization of 70- 
deg. sphere cone, a forebody heatshield thermal protection system of SLA-561, a supersonic disk-gap-band 
parachute, and autonomous terminal descent propulsion. These Viking heritage EDL systems, the thin Martian 
atmosphere, and small scale height of obstacles on the ground limit accessible landing sites to those below - 
1.0km MOLA. So far the southern hemisphere has been largely out of reach. (Approximately 50% of the planet 
surface remains inaccessible with current EDL technologies.) The highest landing to date is MER-Opportunity 
at Meridiani Planum (-1km MOLA). Mars Science Lab (MSL) is attempting to develop an EDL system capable 
of delivering 0.8MT to +2km MOLA. 

The Mars Exploration Program is attempting to develop systems that deliver 0.8 metric tons to the surface 
for MSL and 1.2 MT for Mars Sample Return. In order to prepare for human scale missions, we will need to 
develop EDL systems that can get 30-100 metric tons down to the surface per landing. Viking heritage systems 
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Alt.(km MO LA) 



Vel. (m/s) 


Figure 1: Phase plot for robotic MSL entry with ballistic coefficient (3 = 

m/(C D A) = 100kg/m 2 and lift to drag ratio L/D = 0 18.[1 1] 


are unable to achieve this goal. The problem can be illustrated using an altitude versus velocity phase chart. 
Figure 1 shows the phase chart for the 
robotic MSL mission. Entry into the 
sensible atmosphere occurs at the right 
boundary. The goal is to attain zero ve- 
locity at the target landing altitude at the 
lower left corner. The track ending in a 
dashed line indicates the trajectory with- 
out supersonic parachute deployment 
that hits the ground (0 MOLA) at 200 
m/s. The red triangular region indicates 
the safe zone for supersonic parachute 
deployment - bounded by high and low 
Mach number, high and low dynamic 
pressure, and minimum altitude require- 
ments. The blue and green zones repre- 
sent domains for subsonic parachute de- 
ploy and subsonic retro-propulsion for 
final landing, respectively. Viking her- 
itage EDL requires a trajectory that en- 
ters this triangle, preferably near the 

center of the high Mach boundary, to enable supersonic parachute deployment. The track ending in a solid 
line in Fig. 1 indicates an MSL trajectory with supersonic parachute deployment. 

Figure 2 focuses on the lower left corner of Fig. 1 overlayed with ballistic trajectories as a function of 
ballistic coefficient ( 3 . Increasing ballistic coefficient causes the vehicle to penetrate deeper into the atmo- 
sphere (lower altitude) for a given velocity. If the area of the entry vehicle is constrained by launch ve- 
hicle diameter from Earth then the bal- 
listic coefficient will increase with in- 
creasing payload mass to Mars. The 
message from Fig. 2 is that a ballis- 
tic coefficient of 150 is a limiting value 
for current technology. The altitude is 
too low to safely deploy a supersonic 
parachute above this value of the ballis- 
tic coefficient. Some relief is attained 
with lifting bodies but the fundamental 
limitation remains - heritage EDL sys- 
tems cannot deliver large payloads to 
Mars with given constraints of launch 
vehicle diameter. Two new EDL sys- 
tem elements are currently being inves- 
tigated to overcome this limitation. The 
first element lowers the ballistic coeffi- 
cient by deploying structure to increase 
area. Inflatable ballutes deployed in the 


Alt.(km MOLA) 


Varying /? from 25-200 kg/m 2 
(chute-less trajectories) 



Figure 2: Detail of Figure 1 showing influence of ballistic coefficient on effec- 
tiveness of heritage EDL systems. When (3 is too large altitude is too low at 
Mach 2 to deploy parachute.[11] 


Vel. (m/s) 
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supersonic regime or prior to entry interface are primary examples of this approach. The second element aug- 
ments aerodynamic drag with retro-propulsion engaged in the supersonic and possibly hypersonic regimes. 

2.1.1 Supersonic Retro-Propulsion — SRP 

Korzun et. al. [12] surveyed experimental and computational efforts addressing use of retropropulsion ap- 
plicable to EDL systems. A single centered jet (emanating from the vehicle axis) and peripheral jets were 
compared. They report that for low thrust coefficients (ratio of jet thrust to freestream dynamic pressure times 
area) the peripheral jet system was superior to the center in terms of total axial force coefficient (aerodynamic 
drag plus thrust). The center jet tends to create an aerospike effect that reduces aerodynamic drag. At large 
thrust coefficients the aerodynamic drag is overwhelmed by the thrust and there is no advantage to periph- 
eral jet placement. As computational analyses mature with objective-function optimization the opportunity to 
minimize retropropulsion system mass through tailoring of nozzle design and vectoring should be realizable. 
Indeed, engaging low thrust coefficient retro-propulsion, even in the hypersonic domain, may be useful in re- 
ducing surface heating if a film cooling system can be adapted from propulsion system infrastructure. The 
creative use of CFD for design in this application is currently stymied by an inability to effectively adapt the 
grid in the jet interaction region - even in a supersonic domain. An example of this challenge will be presented 
in the Numerical Results. 


2.1.2 IRDT / Ballutes 


A ballute (balloon-parachute) or Inflatable ReEntry and Descent Technology (IRDT) [13] is an inflatable, aero- 


dynamic drag device for lowering the ballistic coefficient of planetary entry vehicles. Ballutes may be directly 


attached to a vehicle, increasing its cross- 
sectional area upon inflation, or towed be- 
hind the vehicle as a semi-independent de- 
vice that can be quickly cut free when the 
requisite change in velocity is achieved. 
Rohrschneider et. al. [14] survey ballute 
technology for aerocapture. 

An example of a trailing ballute[15, 16] 
is shown in Fig. 3 including the computa- 
tional grid and color contours of pressure 
in the quarter domain symmetry planes. 
(Conditions of the simulation are for a Ti- 
tan Organics Explorer with velocity equal 
to 8.55 km/s and density equal to 1.9 x 
10 -7 kg/m 3 . [17, 18] The gas model in- 
cludes molecular nitrogen and atomic ni- 
trogen in thermochemical nonequilibrium.) 
An inflatable tether (required for compres- 
sive loads encountered during deployment 
prior to atmospheric interface) connects the 
towing spacecraft (lower right corner) to the 
inflated toroid. The bow shock of the space- 
craft is enveloped within the hole of the 



Figure 3: Unstructured grid used in simulation of spacecraft - tether - bal- 
lute system. Colors correspond to pressure levels nondimensionalized by 

PooVl. 
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toroid, interacting with the toroid bow shock. A high pressure zone in the core of the toroid drives flow back 
toward the base of the towing spacecraft. 

Ballute simulation provides challenges on many fronts. Inflatable surfaces deflect and stretch under aerody- 
namic load and relatively high temperature. [19, 20] Unsteady interactions may follow from the surface deflec- 
tions or from a purely fluid dynamic origin (discussed later in more detail) and the ability to accurately capture 
these phenomena are a strong function of the ability to resolve free shear layers and the triple points of shock- 
shock interactions. A wide variation of Knudsen numbers may exist in the same simulation - in the example of 
Fig. 3 the Knudsen numbers vary from 0.026 on the toroid to 0.88 on the tether. This range indicates a variation 
from continuum to near free molecular flow and requires analyses appropriate for these domains in the same 
simulation. [21] The inflation infrastructure provides opportunity to engage film cooling and raise acceptable 
aerothermal loads with innovative material concepts. [22] 


2.1.3 Energetics 

Harvesting and storing energy imparted to the gas during atmospheric entry (regenerative aerobraking) offers 
significant potential to do useful work after landing. Macheret et. al. [23] describe a concept for a regenerative 
aerobraking system at Mars using a magnetohydrodynamic (MHD) system. Though not usually considered an 
integral part of planetary vehicle aerothermodynamics this potential application illustrates that MHD simulation 
should not be ignored in establishing a complete entry vehicle design competency. 


2.2 Exploration of Outer Planets 

A NASA systems study for aerocapture at Neptune reveals significant challenges to the thermal protection 
system design. [24, 25] The study concluded that aerocapture at Neptune with inertial entry velocity of 29 km/s 
could be accomplished with a rigid aeroshell configured as a flattened ellipse-sled with a lift-to-drag (L/D) ratio 
of 0.80 and a ballistic coefficient j3 « 895 kg/m 2 . Laub et. al. [25] note that the large convective and radiative 
heating rates on the nose (10-15 kW/cm 2 ) and on acreage windside (3-6 kW/cm 2 ) using existing materials 
require thicknesses that currently cannot be manufactured as a high quality product. Neither the effects of 
coupled ablation and radiation nor the effects of significant shape change on aerodynamics were considered in 
this high level study but would be required as designs mature - a goal of CA development as discussed below. 


3.0 ALGORITHMS 

In the previous chapter, Problem Definition, we discussed simulation capabilities required for defining the 
aerothermal environments of planetary entry vehicles and opportunities for exploring innovative concepts for 
EDL systems. In this chapter we highlight two algorithmic topics that must be addressed to achieve these 
simulation capabilities. We note that these topics are of special concern, if not unique, to hypersonic flow 
simulation. The first addresses the need to obtain accurate heating simulation on grid systems that can be 
easily adapted to complex flow topologies. The second addresses the infrastructure needs to support fully 
coupled ablation, radiation, and shape change in a CA simulation. Governing equations and physical models 
for atmospheric entry have been presented previously [26] and will not be repeated here. 
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3.1 Heating on Unstructured Grids in Hypersonic Simulations — Multi-Dimensional Recon- 
struction 

The simulation of hypersonic flows with fully unstructured (tetrahedral) grids has severe problems with respect 
to the prediction of stagnation region heating. [27, 28] The problems arise for two reasons. [29] First, good 
shock capturing in the hypersonic regime requires alignment of the grid with the captured shock. Alignment 
here means that all surfaces of a control volume in contact with the shock are either parallel or orthogonal to 
the discontinuity. Any skewness present in a conventional upwind scheme (in which first-order reconstruction 
is a function of only two nodes on opposite sides of a shared face) results in a smearing of the shock, which 
manifests itself as a non-physical distribution of entropy from streamline to streamline crossing the shock. In 
the stagnation region, there is very little dissipation of these entropy gradients as they pass from the shock to 
the boundary-layer edge, and the surface heating is very sensitive to these non-physical gradients. Tetrahedral 
cell topologies, by their very nature, will always have at least one surface that is skewed to the captured shock. 

Second, tetrahedral grids are not well suited to the preservation of symmetry and the biased tetrahedral 
grids used in the simulations of heating herein are particularly ill-suited for this purpose. Asymmetric elements 
and faces skewed to the symmetry surface tend to induce non-physical cross flows. Again, this problem is 
particularly evident in the stagnation region where the physical flow velocities themselves are small and the 
non-physical, grid-induced velocities are of similar magnitude. Surface heating rates are a sensitive indicator of 
these problems. The combination of random jaggedness in the shock capturing process as well as non-physical 
cross flow velocities cross-couple to produce unacceptable surface heating predictions. 

Problems with heating can be overcome if one uses a semi- structured grid (e.g., prisms) across the boundary 
layer and adapts the grid to the shock. Nompelis et. al.[30] show excellent heating results on families of grids 
where a prismatic grid was adapted to the shock and acceptable heating results when an unbiased (random face 
orientation), tetrahedral grid was adapted to the shock. Prismatic elements are relatively easy to generate on 
blunt body geometries and their superior performance with respect to aeroheating predictions have led to their 
use as standard practice in unstructured simulations. There is no better alternative with todays algorithms than 
to use prismatic elements to capture the boundary layer and bow shock. 

If one is doing a hypersonic simulation without any free shear layers or internal shocks then the specialized 
application of prismatic elements at the body and bow shock is perfectly acceptable. If such flow structures 
exist — as they do in almost all interesting problems noted in the previous chapter — then the accuracy of 
any algorithm must be questioned when it ignores special topological grid requirements where the features are 
harder to resolve. Consequently, overcoming these accuracy issues on tetrahedral grids is considered one of the 
highest priority tasks for achieving simulation for innovation. 

A three-dimensional reconstruction algorithm to address these challenges was originally described in ref- 
erence [29]. Some details of the algorithm are included in the Appendix for convenience of the reader. The 
key element of the algorithm is that inviscid flux reconstruction is element-based rather than edge-based so that 
even the first-order reconstruction is based on 4 nodes surrounding the element rather than 2 nodes defining an 
edge. We note that this cylinder challenge problem and a related sphere challenge problem (not included here) 
have been presented in [31]. 

3.1.1 Challenge Problem - Biased, Tetrahedral Grid on Cylinder with Spanwise Resolution 

A challenge problem is identified in which a highly biased, tetrahedral grid is applied to the simulation of a 
hypersonic flow over a 3D cylinder at 5km/s, 0.001kg/m 3 and T ^ 200K in air (M^ « 

17). [28, 32] (In fact, other hypersonic free stream conditions are acceptable. The challenge here is to recover 
the proper spanwise uniformity within 1%.) A structured-grid solution generated with LAURA[33] is used 
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as both a benchmark and to generate initial grids for use in FUN3D.[28] A successful simulation requires a 
constant spanwise simulation of heating, pressure, and shear and symmetric simulation of these quantities on 
the right and left sides of the cylinder. The grid bias accentuates any issues with the algorithm. To date there 
have been no successful simulations of this challenge problem on the given grid using conventional upwind 
formulations. (A recent thesis on a related unstructured grid without bias produced excellent heating symmetry 
using a high-order Discontinuous Galerkin finite element method with a PDE based artificial viscosity. [34]) 
The structured grid from LAURA, 
adapted to the shock and boundary layer, 
is converted from hexahedral elements to 
tetrahedral elements by adding diagonal 
edges consistently from minimum index to 
maximum index corners. A comparison of 
the grids in a plane of nodes perpendicular 
to the cylinder axis is presented in Fig. 4. 

The structured grid has 65 nodes from the 
body to the inflow boundary ahead of the 
bow shock and 61 nodes from the left to 
right outflow boundaries. Node placement 
is identical between the two grids. The 
placement of additional edges in the un- 
structured grid is the only difference. The 
strong biasing of diagonals in the grid is 
an intentional characteristic to expose al- 
gorithm deficiencies that may otherwise be 
averaged out. 


FUN3D 


LAURA 



Figure 4: Structured grid (LAURA, bottom) and biased, unstructured grid 
(FUN3D, top) in plane orthogonal to cylinder surface. 




(a) Structured, quadrilateral surface elements 


(b) Unstructured, triangular surface elements 


Figure 5: Surface mesh on cylinder with ten rows of spanwise cells. 


A key element of this test problem is the addition of ten spanwise cells, shown in Fig. 5, across the cylinder, 
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providing additional degrees of freedom in the simulation to allow asymmetries to develop. The span wise 
degrees of freedom enable non-physical cross-flow velocities to develop and irregular shock capture in the 
spanwise direction that combine to corrupt the predicted heating distribution. Earlier tests had shown that 
single spanwise cell simulations were not severely compromised. 



(a) Heating (blue squares) and shear (red squares) distribution 
for each of ten rows of spanwise cells. 
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Circumferential 


(b) Surface heating contours. 


Figure 6: Simulation results using conventional reconstruction on the tetrahedral grid with ten spanwise cells. 



(a) Values of heating - blue circles, shear - red triangles, and 
pressure - black squares at nodes on surface. Figure shows results 
for all 1 1 surface nodes at a given x location. 


q, W/cm 2 : 10 1520 2530 35 40 45 50 5560 



Circumferential 


(b) Heating contours on surface showing significant improve- 
ment in spanwise symmetry. Some edge effects at spanwise 
boundaries still persist. 


Figure 7: Simulation results using Option 2 reconstruction on the tetrahedral grid with ten spanwise cells. 


A large spanwise variation in heating and shear for this simulation using conventional reconstruction is 
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evident in Fig. 6(a) with ±30% about the mean — an unacceptable result for aerothermodynamic analyses. 
Recall that the symbols at every x location should overplot and the peak value of shear on the right should equal 
the peak value of shear on the left. The LAURA result in this case is 52 W/cm 2 which is in good agreement 
with the spanwise mean of the heating at the stagnation point. The surface heating contours in Fig. 6(b) indicate 
that a slight drift in the velocity — probably induced by grid bias - produces a higher heating toward the front 
(y = 0) plane. 

Results for the challenge problem conditions are repeated in Fig. 7 using multi-dimensional reconstruction. 
The new results (using Option 2 as defined in the Appendix) are presented as symbols and a benchmark, 
grid-converged, structured grid solution (LAURA) for the same case are presented as solid lines. The heating 
contours show excellent spanwise uniformity in Fig. 7(b) in the interior. Some end effects are evident at the 
top and bottom (spanwise) boundaries. The heating and shear distributions are in good agreement with the 
benchmarks. 

The element-based reconstruction algorithm requires three reconstructions per element whereas the base- 
line ID algorithm requires one reconstruction per edge. The cylinder challenge problem with 10 spanwise 
cells required 14.8 seconds per relaxation step as compared to 6.8 seconds per step for the baseline ID algo- 
rithm. Some improvement is expected when the inviscid and viscous terms are computed in a single loop over 
elements. This refactoring will eliminate some redundancies occurring in both the inviscid and viscous formu- 
lations. Additional simulations using multi-dimensional reconstruction will be presented in the next chapter on 
more challenging problems. 

3.2 Infrastructure — Free Energy Minimization - (FEM) 

The Gibbs Free Energy minimization (FEM) option is implemented in LAURA with a goal of making simple 
adjustments to the baseline thermochemical nonequilibrium flow solver to enable a thermochemical equilibrium 
solver for generic mixtures of gases with spatially varying elemental mass fractions. The capability is required 
to enable interfaces with ablation routines based on equilibrium gas chemistry. The capability also enables 
simulations that are so deep into the equilibrium regime that the magnitude of species chemical source terms 
overwhelm the magnitude of convection and diffusion. (Roundoff error in the evaluation of the source terms is 
significant relative to the magnitude of convective and diffusive terms.) Consequently, the sum of the species 
continuity equations does not telescope to a meaningful mixture continuity equation because of roundoff error. 
The modifications are originally described in [35]. In essence, the solution of N s species continuity equations 
with the solution of N e elemental continuity equations and ( N s — N e ) algebraic, partial equilibrium equations. 


3.2.1 Elemental Continuity Equations 


The transformation from species continuity to elemental continuity involves the identification of the fraction 
of element i in each species j. The ~ overset symbol denotes an elemental property as opposed to a species 
property. The transformation matrix F involves N e rows, one for each element, and N s columns, one for each 
species. The fraction is computed by multiplying the elemental weight, by the number of atoms of that 
element in the species, and dividing by the species weight Mj. 


F — 

r h3 ~ 


&i,j 

Mo 


Consider the one-dimensional species conservation equations set 

dq df 

( — w 

dt dx 


CD 

( 2 ) 
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where f includes convection and diffusion. If one multiplies both sides of Eq. 2 by F one obtains the elemental 
conservation equations 


dq df 

dt + dx 
<9q df 
dt + dx 


Fw 

0 


(3) 

(4) 


Note that Fw = 0. Chemical reactions cannot create or destroy an element in the control volume. Conse- 
quently, there is no need to add the source term to the species continuity equations prior to assembling the 
elemental continuity equations. Besides saving operation counts, it is important to apply the transformation 
matrix before adding the source term because the critical components of continuity through convection and 
diffusion may be of the order of roundoff error in the assembly of the source term w. The assembly of the 
Jacobian of the elemental conservation laws (Eq. 4) employs the same matrix transform and retains the use of 
species densities (rather than elemental densities) as the fundamental dependent variables. Consequently, the 
remaining non-base species must be solved using partial equilibrium relations in order to close the equation set 
for species densities. 


3.2.2 Partial Equilibrium Relations 


The derivation of partial equilibrium relations through the minimization of free energy is available in several 
texts. [36] Consider the partial equilibrium relation for non-base species j which may be expressed as 

reactants products 

Y. < = > y Pn,jN n (5) 

m n 


where one of the products N n corresponds to non-base species j, the reactants correspond to constituent base 
species, and one additional species may be required to enforce a chemical balance. The partial equilibrium 
relation involving non-base species j may then be derived as follows. 
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where R = p^R/ 101, 325 to include conversion from pressure in atmospheres to the metric system and to 
non-dimensionalize density by p^. The non-base species continuity equations are thus replaced by a linear 
equation in In pp with coefficients defined 


{ 0 if species k is not part of reaction j 

—/3 n ,j if species k — n is a product in reaction j 
/3 m j if species k = m is a reactant in reaction j 


( 10 ) 
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The temperature dependent equilibrium constant Kj ( T ) is curve fit using 

Kj (T) = A A + Bj + Cj In Z + DjZ + EjZ 2 (1 1) 

Zj 

where Z — 10, 000/T and the constants Aj - Ej are defined by fitting to the left-hand- side of Eq. 8 (defining 
Kj (T) in the left hand side of Eq. 9) at T = 1000, 2000, 4000, 8000, 16, 000. The enthalpy H n and entropy 
S n for species n is obtained from curve fits. [37] 

The coupled solution of Eq. 9 with the other conservation laws can be a challenge, particularly when it 
is early in the simulation and the solution is far from a converged state. It has been found that ignoring the 
Jacobian of Eq. 9 with respect to temperature and by multiplying the Jacobian of Eq. 9 with respect to species 
densities by a safety factor of 15 has enabled solution of the FEM coupled equation set for all tests conducted 
for this paper. The factor 15 was determined through numerical experiment. 

3.3 Infrastructure — Ablation Boundary Conditions 

The ablation boundary condition derived here assumes there will be coupling with an ablation module such as 
FIAT[38] in conjunction with an equilibrium chemistry module. It is assumed that the ablation mass flow rate 
rh = rhpyroiysis + ™char, the elemental composition of ablation gases c^ a &/, and the temperature of the ablating 
surface are supplied by the ablation module or can be user specified. The ablation module itself requires the 
aerothermodynamic environments (pressure, radiative and convective heating, and possibly shear) over time to 
compute material response. Consequently, a boundary condition must be applied to the flow solver that closes 
the equation sets for either nonequilibrium external flow using species continuity equations or equilibrium 
external flow using elemental continuity equations. As noted in the previous section, the species densities are 
retained as the fundamental dependent variable data set regardless of the assumption of equilibrium state of the 
external flow. 

It is assumed that the gas chemistry at and within the ablating heatshield is in an equilibrium state. A steady, 
one-dimensional, elemental continuity equation (Eq. 12 ) describes the flow of an equilibrium gas across the 
boundary that accounts for convection and diffusion of element i. Here again, the ~ overset symbol denotes an 
elemental property as opposed to a species property. 


dpjv _ dpjVj 
dy dy 


( 12 ) 


Because we wish to work directly with species densities the transformation matrix introduced in the previous 
section is applied in Eq. 13 below. 


dpjv 
,J dy 



(13) 


where summation over species j is implied. 

Integrate Eq. 13 over y from some virtual point sufficiently deep in the heatshield where it is assumed that 
the net elemental mass flux is given by c^/ra. This assumption ignores the possibility that ablating mass flux 
sources are distributed through the depth of the char. Short of simulating flow through the char, it is not clear 
that this modeling assumption can be improved. The resulting relation (Eq. 14) provides a balance for each 
element in which the convected elemental mass flux out of the boundary minus the diffused elemental mass 
flux into the boundary must equal the ablated elemental mass flux specified for the material. 
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Fijpjv ~ Ci iM rh = FijpjVj (14) 

d Xj 

Fijpjv - c^abim = F l , J pDj— (15) 

Note that an effective binary diffusion coefficient, Dj , is applied to model diffusive flux. A mass corrected 
diffusive flux[39] is also applied to guarantee that the sum of diffusive fluxes over all species (or over all 
elements) goes to zero. Simulations show some elemental separation near the wall for the no blowing case. 
Elemental separation in the no-blowing case is less when modeling multi-component diffusion using the Stefan 
Maxwell equations. [39] A verification check shows that no elemental separation is computed in the no-blowing 
case when a constant Schmidt number is applied. 

The discretization of Eq. 15 is shown in Eq. 16. It is necessary to solve for all dependent variables in the 
ghost cell 0. 


\pjV)l + (pjv) o 


h3 


Ci,ablFl — Fij 


(pDj)i + (pDj) 0 


Xj, l ~ Xj , o 

Ay 


(16) 


The solution in the ghost cell requires simultaneous relaxation of the partial equilibrium relations given by 
Eq. 9. The Jacobian formed using In pj as the dependent variable provides a more robust convergence of this 
boundary condition. A relaxation factor is applied to each row of the Jacobian according to the magnitude of the 
current residual for that row equal to n/(d In [pj] 2 +n) where n is the local iteration count on relaxation steps for 
this boundary condition. Note that in the limit of zero blowing this boundary condition specifies an equilibrium 
catalytic wall - that is, the species are in chemical equilibrium as determined by the wall temperature and 
pressure and elemental mass fraction gradient equal to zero. 

\ [(Hi + (Ho] = ™ ( 1? ) 


Pi + pivj = Po + poVo (18) 

The equation set is closed with a discretization of net mass flux (Eq. 17) and normal momentum balance 
(Eq. 18). Eq. 17 is used to eliminate vq f rom all equations. Then Eq. 18 replaces one of the elemental continuity 
equations so that the ghost cell pressure po is consistent with the specified temperature Tq = 2 T w — T\ and the 
computed densities p 3 p. 

At present, ablation rates and ablation gas elemental mass fractions are explicitly set by the user. The 
boundary condition is sufficiently robust to enable a coupled simulation assuming equilibrium gas chemistry at 
Voo =15 km/s and a dimensionless blowing rate equal to 0.2. 

On the first step of initialization, all cells (including ghost cells) are set to inflow conditions. On the second 
step of initialization, ghost cell velocities are reset to zero and wall temperatures are set to a user specified 
value. Simulations involving ablation and/or radiation are started from a converged (or nearly converged) non- 
radiating and non-ablating solution. The line-implicit formulation for relaxation of interior cells requires a 
Jacobian of the boundary condition. It is assumed for this purpose that species mass fractions and temperature 
in the ghost cells are fixed. (The point-implicit formulation for relaxation of interior points, by definition, 
requires no Jacobian information of neighbor cells.) 

During interior cell updates the surface convective heating, shear, and diffusion flux are calculated from 
current values in the ghost cells. Ghost cell values, in turn, are updated at a user defined frequency, N exc h an ge- 
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(Ghost cells are also used to hold information from neighbor blocks and the variable name N exc h ang e derives 
from the frequency in which neighbor block information is exchanged.) Ghost cell values of species density 
are fully defined by interior point values - there is no time dependent update or damping using any fraction of 
a previous value. Damping is applied to the update of a radiative equilibrium surface temperature in the case 
of no blowing. The radiative equilibrium wall temperature boundary condition is then defined with T w new — 
e\q/ (ecr )] 1 / 4 + (1 — e)T w old . In this case, the ghost cell temperature is computed as Tq = 2 T w — T\. 

N exchange is set equal to 2 in all cases presented herein. The temperature damping factor e is varied between 
0.001 and 0.01. In more recent applications (some involving a surface energy balance not covered here) the 
combination N exc h an ge — 20 and e — 0.1 have provided more robust convergence. 

3.4 Infrastructure — Radiation Coupling 

Radiative energy transport is computed using the tangent slab approximation with the HARA code. [40, 41]. 
The flow solver provides species densities, temperatures, and distances along a ray orthogonal to the surface 
as inputs to HARA. HARA returns divergence of the radiative flux and radiative flux to the wall for use in 
LAURA. Calls are made to HARA on the order of hundreds of relaxation steps, N ra d , holding the radiative 
quantities constant over this period (loosely coupled). Options are provided to interpolate radiative quantities 
between selected rays to reduce overall time in the radiation solver. The HARA radiation model treats atomic 
lines and continuum of N, O, C, and H. Molecular bands are treated for N 2 , N^, 0 2 , NO, C 2 , C 3 , CN, CO, C 2 H 
with the smeared-rotational band (SRB) method for molecular band radiation. Non-Boltzmann electronic state 
models are applied for atoms and molecules. Typical computational time for a single line-of-sight is 20 to 50 
seconds - short enough to enable fully coupled simulations with codes like LAURA. 

A parameter f ra( i may be specified to introduce the divergence of the radiative flux more slowly into the 
flow solver. 


V • <lrad = fradV ' Qr ad (HARA) + (1 - f rad )V • q? ad (19) 

If the radiation introduces strong cooling into the shock layer the shock standoff distance decreases signif- 
icantly. As the captured shock moves closer to the body, computational cells previously in the hot shock layer 
are overtaken by the pre-shock domain. However, these cells retain the cooling term V • q ra d,i,j,k from the 
previous call to the radiation solver. The LAURA code enforces a floor on temperature for every point in the 
flowfield that allows the simulation to survive such transients under most circumstances. However, cases have 
been observed where introduction of the full update to radiation has caused stability problems. Best practice so 
far has been to set f ra g to 0.4 on the first call from an uncoupled solution and then to increase its value to 1 on 
subsequent returns from the radiation solver. 

4.0 SIMULATIONS FOR VALIDATION 

4.1 Validation Challenges 

Let us focus on the unique challenges of code validation for simulation of hypersonic flow for planetary entry. 
The key issue in the hypersonic flight domain is to adequately model the mechanisms through which a massive 
reservoir of kinetic energy cascades across diverse chemical and internal modes. More precisely, the kinetic 
energy of an entry vehicle cascades into chemical, translational, and internal modes of the quiescent atmospheric 
gas through action of the vehicle’s bow shock. The transfer of energy through space by laminar and turbulent 
transport, convection and radiation are related critical elements in a standard simulation. Add to this mix the 
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physics of gas-surface interactions (catalysis), material response (ablation), jet interactions (retro-propulsion), 
and/or elasto-dynamic boundary conditions (inflatable / deployable structure) and the challenges to planetary 
entry validation are daunting. 

The path forward must combine measured data from ground-based tests and flight experiments as well as 
computed data from first-principles physics simulations on canonical problems (computational chemistry [42], 
DNS[43]). Each has pros and cons. We confine our attention here to measured data. 

Ground-based tests provide the best opportunity to measure boundary conditions and probe the flow field 
with sufficient detail to validate physical models. However, the facilities are not generally quiet and often the 
flows are not perfectly uniform nor steady. As enthalpies are increased, flow quality and test times become 
even more problematic. Recent efforts by Candler et al. [44] underscore the proposition that future code 
validation tests will require more complete simulation of the facility flow in addition to simulation around 
models. Simulation of facility flow with requisite validation checks should become a standard component 
of CA analysis practice. [45] One thus acknowledges that ground-based facilities often “bruise” the gas (flow 
uniformity, chemical and thermal state) but in fact these perturbations provide additional opportunity to validate 
the suite of physical models used in CA simulations. 

As noted by Wright el. al. [46] “Flight testing, while certainly capable of fully exercising the TPS design, 
is prohibitively expensive in most cases. In general flight testing should be reserved for final model and system- 
level validation, once we have good physics-based models of the expected environment and resulting material 
performance. Ideally, the purpose of a dedicated flight experiment should not be to learn anything new about the 
flight environment encountered, but rather to validate that the models employed are adequate.” However, valu- 
able engineering information can be obtained at relatively lower cost by including engineering instrumentation 
on science missions. While there is little benefit to the mission that flies it, the overall exploration program ben- 
efits from the additional knowledge gained about the flight environment. This new knowledge may then benefit 
future missions through reduction of TPS mass requirements, for example, thus enabling increased payload 
mass. The Mars Science Laboratory Entry, Descent, and Landing Instrumentation (MEDLI) Project is a prime 
example of piggybacking engineering instrumentation on a science mission - Mars Science Laboratory. [47] 
The flight science objectives of MEDLI directly address the largest uncertainties in the ability to design and 
validate a robust Mars entry system, including aerothermal, aerodynamic and atmosphere models, and thermal 
protection system (TPS) design. [47] 

The remainder of this chapter provides example of code assessment and validation using the structured flow 
solver LAURA and the unstructured flow solver FUN3D in problems of particular interest to the CA community. 
Both ground-based and flight data are simulated. As CA practice matures the simulation of each of these cases 
should become routinely accessible with published uncertainties of predicted aerothermal environments. 

4.2 Retro-Propulsion 

Simulations are executed to assess the ability of CFD to predict the flow field associated with a jet firing from 
a blunt body into on oncoming supersonic flow. Of particular interest is the ability to achieve a grid converged 
solution in a flow topology where the location of critical domains requiring fine resolution are not known apriori. 
The sensitivity of flow topology and aerothermal loads as a function of turbulence modeling are also quantified 
in these studies. The test case involves a Mach 2 oncoming flow. ( Voo = 527m/s, = 3.55 x 10 -2 kg/m 3 ) 

[48]. Experimental data indicates this condition results in a steady interaction - removing one significant level 
of complexity from the simulation. 

The simulations in Fig. 8 used the unstructured grid solver FUN3D. The model is a spherically capped 
cone with a cylindrical afrebody. The nose of the vehicle and exit plane of the nozzle is seen in this figure 
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(a) Turbulent shear layer - SST model 

-3 



(b) Turbulent shear layer - Spalart - Almaras model 

o 1 



(c) Laminar flow 

Figure 8: Mach contours in simulation of retro-propulsion in a Mach 2.5 oncoming flow. 


at x = 0. The initial grid generation process was found to be critical in this application. Grid adaptation 
driven by local gradient information failed to produce grid converged solutions. Rather, the grid converged on 
intermediate flow topologies and then locked in on these states - failing to evolve through a coarser grid that 
bounds the high gradient regions of the intermediate states. Ultimately, the expert user was required to throw 
a relatively fine initial grid across the entire interaction domain with judicious use of source terms in the grid 
generation procedure before a grid converged solution with qualitative agreement with experimental data could 
be obtained. (This situation is repeated in the next test problem - the sharp double cone. A fine initial grid 
is required to enable the separation zone to progress upstream as far as it needs.) Approximately 4.3 million 
nodes (25 million tetrahedra) were used to resolve one quarter of the flow domain surrounding the model (zero 
degrees to ninety degrees around the symmetry axis). The requirement here for an expert in the loop (in this case 
an expert in training!) to recognize that the grid adaptation strategy was converging on a wrong solution is a 
show stopper to simulation for innovation. How can one explore innovative concepts for thruster configurations 
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when the simulation for any one configuration requires so much oversight? The application of adjoint-based 
adaptation[49-51] is expected to overcome this problem of convergence to the wrong answer and thus is a high 
priority task on the path to simulation for innovation. 

With the requisite initial grid generation strategy established, it was possible to study effects of turbulence 
modeling on the interaction. These scoping results using the SST model [52] (Fig. 8(a)), the Spalart Almaras 
model [53] (Fig. 8(b)), and laminar flow (Fig. 8(c)) show significant variation of flow topology. Indeed, the 
variation is not unlike that seen with the different initial grids discussed above. Results using the SST model 
are in closest qualitative agreement with the experimental data. Significance of agreement is only qualitative 
because only Schlieren results were available and some inflow conditions were not defined. A new experimental 
validation program has been initiated and supporting simulations are being used to assess possibility of flow 
blockage. 

4.3 Sharp, Double Cone 

A series of experiments to measure pressure and heating for code 
validation involving hypersonic, laminar, separated flows over dou- 
ble cones was conducted at the Calspan-University at Buffalo Re- 
search Center (CUBRC) in the Large Energy National Shock (LENS) 
tunnel. [54] The configuration (Fig. 9) yields a flow topology with a 
large separation zone across the intersection of the two cones. An 
overview of the flow topology is evident in the temperature and pres- 
sure contours of Fig. 10. The shock off the forward part of the sepa- 
ration intersects the bow shock over the second cone. A transmitted 
shock from this intersection reflects off the wall to a contact surface. 

The extent of separation and the location of pressure maxima down- 
stream of the reattachment point are sensitive to gas chemistry. These 
features are also a sensitive function of grid convergence. This sen- 
sitivity to physical models and grid convergence makes simulation of 
this flow an excellent code validation test. The experiments were the 
focus of a blind code validation study. [54] A single case, Run 28, is featured here with inflow conditions 
Vqo = 2.664 km/s, = 6.55 10 -4 kg/m 3 , T \ ^ = 185.6 K, and wall temperature T w = 293.3 K. These 
conditions correspond to = 9.59 and Re ^ = 144, 010 /m in Nitrogen (to minimize effects of chemical 
reactions). The conditions are the same as originally presented to the study participants. A subsequent evalua- 
tion of test conditions was spurred by a consistent disagreement (20%) between experimental data for heating 
on the front cone prior to separation and multiple, grid converged simulations. [55] This study found that: (1) 
vibrational freezing of flow exiting the facility nozzle; (2) vibrational energy accommodation at the surface; 
and (3) small levels of flow non-uniformity must be included in the simulations to improve the heating predic- 
tion. These corrections to boundary conditions are not included in these results. The intent here is to show 
that the simulations on a tetrahedral grid system are consistent with the earlier, grid converged, structured-grid 
simulations. 

The computed pressure coefficients and Stanton number distributions from three simulations are compared 
to experimental data in Fig. 11. All three simulations use the same set of nodes and Roe / STVD based inviscid 
flux formulation. The LAURA solution was generated for the blind validation study (solid blue line)[56] and 
uses a cell-centered, structured-grid formulation. The FUN3D simulation on a hexahedral grid (small black 
circles) uses a conventional, one-dimensional, edge-based inviscid flux formulation. The FUN3D simulation 



Figure 9: Dimensions (cm) of sharp, double 
cone. 
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(a) Temperature contours. 



(b) Detailed view of pressure contours and streamlines in separa- 
tion region. 


Figure 10: Flowfield over sharp, double cone. 


C p , FUN3D, Hex Grid 

C p , FUN3D, Tet Grid (3D Recon.) 


4 - C p , LAURA 



x/L 



x/L 


(a) Pressure coefficient. (b) Stanton number.. 

Figure 11 : Comparison of experimental data with computation for the sharp, double cone, Run 28. 


on the tetrahedral grid (small green triangles) uses the new three-dimensional, element-based inviscid flux 
formulation. The 3D reconstruction on tets indicates a slight increase in the extent of separation. The FUN3D 
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hex simulation falls almost exactly under the original LAURA simulation. The primary differences between 
the two hex simulations is that LAURA is cell-centered and uses a directionally dependent entropy fix while 
the ID reconstruction used for hex grids on node centered FUN3D uses a directionally independent entropy fix. 
The peak pressure appears beneath the transmitted shock at x/L « 1.3. All three simulations begin the rise to 
peak pressure slightly after the experimental data (in red squares). This delay is indicative of the slightly larger 
extent of separation with the original boundary conditions. The reflected waves in all three simulations closely 
follow each other downstream of the peak pressure and heating. This result demonstrates the consistency of the 
3D reconstruction on tetrahedral grids with both structured and unstructured grid simulations using hexahedral 
grids for a complex flow in which principal directions vary significantly across the domain. 

All three simulations employed local-time- stepping (constant CFL). When starting from a uniform flow 
this process produces a fully attached flow (no separation) with shocks set close to the surface in the early 
stages of the simulation. As the relaxation process continues, the shocks move out from the body and the 
primary separation bubble forms at the intersection of the two cones and slowly grows. In both the LAURA 
and FUN3D hexahedral grid simulations, the separation bubble continues to grow until it reaches its converged 
state. In the 3D reconstruction on tetrahedral grids the growth of the separation bubble followed much the same 
history. However, when the leading edge of the separation bubble reached x/L ^ 0.7 the bubble shape became 
unsteady, giving way to a pulsing that could not be quenched with any variation of CFL or entropy fix. A steady 
solution could only be recovered by engaging a time accurate simulation with local sub-iterations employed 
between time steps. The use of time accuracy had been required by LAURA for a related flow condition at a 
higher Reynolds number (Run 24) of the original study. The 3D reconstruction appears to be more sensitive to 
this phenomena. 

Expert oversight requirements for this problem are somewhat counter-intuitive. The LAURA code has two 
grid adaptation features that were tested on the initial simulation. The first feature redistributes points across 
the boundary layer and shock layer according to a user specified cell Reynolds number criteria and percentage 
of points requested in the boundary layer. This first feature worked flawlessly on this simulation. The second 
feature has a bow shock sensing algorithm that adjusts the length of the computational line of points (structured 
grid) such that the inflow boundary is aligned with the captured bow shock. This feature actually caused more 
problems with the evolution of the simulation. The shock transmitted off the leading edge of the separation 
bubble caused significant readjustment of the outer boundary. Each grid update induced changes in the separa- 
tion point. The feedback process — change in separation point to change in outer boundary and grid position to 
change in separation point ... — sustained large, pulsing motion across the cone. Ultimately, the best solution 
was to apply uniform grid refinement across the entire domain until a grid converged solution was obtained and 
forego any attempt to save grid resources with shock alignment. These observations suggest this problem is an 
excellent test bed for grid adaptation algorithms ultimately required in simulation for innovation. The uniform 
refinement approach required a grid of approximately 512 x 1000 cells to demonstrate grid convergence for the 
256 x 500 solution. An intriguing challenge is to execute a grid converged solution with robust adaptation using 
an order of magnitude fewer points. Here again, application of adjoint-based adaptation needs to be proven on 
such a challenge problem. 


4.4 STS - 2 

CFD simulations of heating over the shuttle orbiter (STS - 2), including effects of finite surface catalysis, 
have been reported for the multi-block, structured grid solver LAURA. [57] Some of these solutions were 
redone with a refined grid and updated physical models in support of the Columbia accident investigation and 
Return to Flight activities. [58] The original simulations used a grid with approximately 1 million cells and a 
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7-species air model with an early model for surface catalysis. The more recent simulations on a refined grid 
with approximately 8 million cells used a more current catalysis model[59] for reaction-cured glass (RCG) and 
a 5-species air model (NO + and e - were ignored to save computational time with focus on relative effects 
due to damage). Two new simulations are run with FUN3D for a Mach 24 case[57]: = 6.920 km/s, 

Poo = 5.750 10“ 5 kg/m 3 , ^ = 202 K , a = 39.4 deg with wall temperature defined assuming radiative 
equilibrium with a surface emissivity of 0.89. The structured grid was converted to FUN3D format and a 
simulation was executed using edge-based reconstruction on hexahedral elements. The hexahedral grid was 
converted into a pure tetrahedral system with identical nodes and a second FUN3D simulation was executed — 
this time using element-based reconstruction on the biased, tetrahedral-grid system. 

Two cuts along the windside of the shuttle orbiter are illustrated in 
Fig. 12 over a flooded contour plot of surface heating. Heating levels 
from the original LAURA simulation from 1993, an updated LAURA[60] 
simulation, and the FUN3D simulation on a hexahedral grid and a tetra- 
hedral grid are compared in Fig. 13 along these cut lines. The updated 
LAURA solution for heating is in excellent agreement with the FUN3D 
simulation on a hexahedral grid. These two simulations use identical 
modules for thermodynamics, transport properties, chemical kinetics, 
and surface catalysis. The heating distribution with these new simula- 
tions are generally in good agreement with the experimental data. These 
simulations did not include a deflected body flap. Though not shown here, comparisons along other cuts across 



Figure 12: Cut lines over windside of 
shuttle orbiter showing surface heating 
contours. 


STS-2 Mach 24.3 STS-2 Mach 24.3 




(a) y=0 in. 


(b) y=185. in. 


Figure 13: Comparison of heating on windside shuttle test case among original LAURA (black line), updated LAURA version 5 
(green line), FUN3D on hex grid(red line), FUN3D on tet grid with 3D reconstruction (blue line), and experimental data (black error 
bars). 


the body (windside, leeside, circumferential) show equivalent trends. The older LAURA results over-predict 
both the flight data and the newer simulations along these windside cuts. It is thought that the updated model 
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for surface catalysis is the primary reason for this difference though a systematic study has not been executed 
to check this supposition. The main point of this comparison is to demonstrate that the node based FUN3D 
provides equivalent results to the cell based LAURA when identical grids and physical models are employed. 




Figure 14: Comparison of surface contours on windside of shuttle. Upper surface to symmetry plane was computed using a 
hexahedral grid system. The lower surface to the symmetry plane was computed using biased tetrahedral grid system. 


The next step is to demonstrate that the FUN3D sim- 
ulation on tetrahedral grids with 3D reconstruction pro- 
vides equivalent surface heating as the hexahedral grid 
simulation. This step has been obstructed by a sensitiv- 
ity to the formulation of the outflow boundary condition 
off the trailing edge of the leeside of the wing and ver- 
tical tail. The boundary condition has been revised to a 
vacuum type boundary (enforce low pressure on exit) if 
flow reversal is detected through the outflow boundary. 
The intent is to nudge a reverse flow to exit the outflow 
boundary. The simulation runs with this fix though some 
irregularities in the heating near the outflow boundary on 
the leeside wingtip are evident. 

Of greater concern is the appearance of a narrow 
band of low heating along the plane of symmetry. In 
Fig. 14(a) we show heating contours on the lower sur- 
face of the shuttle orbiter. The upper half of the figure 
was simulated with FUN3D using a hexahedral grid sys- 
tem. The lower half of the figure was simulated indepen- 
dently with a biased, tetrahedral grid system using iden- 



Figure 15: Comparison of surface heating contours on lee- 
side of shuttle. Lower surface to symmetry plane was com- 
puted using a hexahedral grid system. The upper surface to 
the symmetry plane was computed using biased tetrahedral 
grid system. 
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tical nodes. Ideally, the bottom half should reflect the top half. In this case, heating contours near the symmetry 
plane on the lower, tetrahedral grid side abruptly run up toward the nose indicating a cooler region near the 
symmetry plane. Further away from the symmetry plane the heating contours on top and bottom appear to be 
in reasonably agreement as confirmed in viewing the blue line plots in Fig. 13(b). Note that pressure contours 
in Fig. 14(b) show reasonable reflective properties across the symmetry plane - indicating again that heating 
and not pressure should be used as a metric of solution quality in hypersonic simulations. Heating distibutions 
on the leeside (Fig. 15), while not perfectly reflective, do not reveal any inherent problems in execution using 
tetrahedral grids. Note the heating irregularities at the leeside wingtip reported earlier. 

The cause of the cool streak down the windside centerline in the tetrahedral simulation is unknown. A 
suspected cause - an asymmetry in the way the eigenvalue limiting is engaged as a function of principal direction 
(see Appendix) - has been investigated by making radical changes in the limiting process. Radical changes have 
also been made in the definition of principal direction. In all cases, the quality of the surface contour pattern 
is not significantly altered. Future tests will involve AUSM formulations in the context of multi-dimensional 
reconstruction. 

4.5 Ground Based Radiation Data (EAST at NASA Ames) 

Johnston [40] presented a comparison between EAST shock tube radiation measurements [61] and the HARA 
radiation model [62, 63]. Equilibrium and nonequilibrium radiation measurements were studied for conditions 
relevant to lunar-return shock-layers; specifically shock velocities ranging from 9 to 1 1 km/s at initial pressures 
of 0.1 and 0.3 Torr. The simulated shock-tube flow was assumed one-dimensional. The shock-tube flowfield 
was simulated using the LAURA Navier-Stokes solver to predict the flow around a 5 m sphere at the specified 
shock-tube conditions. This multi-dimensional simulation provides a simple way in the context of LAURA to 
generate a standing shock and focus on the quasi-one-dimensional relaxation immediately behind the normal 
bow shock. A two-temperature thermochemical nonequilibrium 11-species (N, N + , NO, NO + , N 2 , N 2 + , O, 
0\ O 2 , 02 + , and e") air model [64, 65] was applied in these computations, and the detailed nonequilibrium 
radiation prediction was obtained from the HARA code in an uncoupled manner (meaning radiative coolong 
was not accounted for). Similar analyses have been performed using other flowfield and radiation solvers by 
Bose et al. [66] and Panesi et al. [67], among others. 

A significant fraction of the radiative heating at typical peak-heating lunar-return conditions is emitted 
from shock-layer regions in thermochemical equilibrium. According to flowfield predictions, the 0.3 Torr 
EAST cases are in thermochemical equilibrium throughout most of the usable test time. Thus, they provide a 
valuable means for validating the present radiation model in the important equilibrium regime. Comparisons 
between experiment and predictions in this regime are also valuable because the problems of nonequilibrium 
chemistry and non-Boltzmann state-populations do not complicate the theoretical modeling, therefore allowing 
the comparisons to focus on the equilibrium radiation properties and not the nonequilibrium rates. 

A comparison between the measured and predicted intensity spectrum at 3.2 cm behind the shock is shown 
in Fig. 16(a) for the 10.34 km/s case at 0.3 Torr. The measured wavelength range for this case is 700 - 900 
nm, with the dominant spectral features being the atomic line multiplets. The significance of these various 
spectral features is indicated in this figure by the cumulative wavelength-integrated ( J c ) curve, which has been 
multiplied by 100 for scaling purposes. Note that a constant value of 0.8 W/cm 3 /sr/^ was subtracted from the 
measurements for this case to force the continuum radiation to agree, as discussed in [40]. This subtracted 
value, which result is a total J c xl00 of 14.4 W/cm 3 /sr, is interpreted as an artifact of calibration or an untreated 
radiation mechanism from contamination species. The J c predictions at shock- velocities 1.5% below the spec- 
ified value are shown in Fig. 16(a) to indicate the sensitivity to the reported uncertainty in the measured shock 
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(a) Intensity as function of wavelength and cumulative intensity (b) Integrated intensity(800-830nm) as function of distance behind 
3.2 cm behind shock front. shock. 


Figure 16: Comparison of measured and computed intensities (700 - 900 nm) behind shock front in EAST facility at 10.34 km/s 
and 0.3 torr. Dashed lines show predictions if specified velocity is reduced by 1.5%. 




(a) Intensity as function of wavelength and cumulative intensity (b) Integrated intensity (850-8 80 nm) as function of distance behind 
3.1 cm behind shock front. shock. 

Figure 17: Comparison of measured and computed intensities (700 - 900 nm) behind shock front in EAST facility at 9.165 km/s 
and 0.1 torr. Dashed lines show predictions if specified velocity is increased by 1.5%. 


velocity. The agreement between the measured and predicted J c values is within 10% throughout this wave- 
length range. Similar agreement is seen in Fig. 16(b), which shows J c values from a sample wavelength range 
as a function of distance behind the shock. The influence of the 1.5% shock- velocity uncertainty is shown in 
these figures along with the prediction obtained by assuming a Boltzmann distribution of electronic states. The 
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Boltzmann and non-Boltzmann predictions are seen to converge in the regions of thermochemical equilibrium, 
as required by theory. 

The prediction of nonequilibrium radiation is dependent upon the modeling of the nonequilibrium chemistry 
and the non-Boltzmann population of the radiating states, in addition to the oscillator strengths and line-widths 
required for equilibrium radiation predictions. The nonequilibrium chemistry and non-Boltzmann models, 
which are assumed to be uncoupled, have significantly larger uncertainties than the equilibrium radiation prop- 
erties. These large uncertainties are due mainly to a lack of experimental measurements at the high-temperature 
conditions of interest in a peak-radiative heating shock-layer. The EAST measurements therefore provide a 
valuable set of data for assessing the non-Boltzmann rate-model assembled by Johnston [63] and applied in 
HARA. These measurements may also be used to validate the quasi-steady state assumption [68] applied in 
the non-Boltzmann models solution procedure, which allows the non-Boltzmann model to remain uncoupled 
from the flowfield model, as opposed to more complicated coupled approaches[69, 70]. Furthermore, be- 
cause the electron temperature and flowfield number densities obtained from the LAURA flowfield solver drive 
the non-Boltzmann model, these nonequilibrium data provide an opportunity to validate the thermochemical 
nonequilibrium models applied by LAURA. 

The study of the 0. 1 Torr data will be similar to the study of the 0.3 Torr data presented above. However, the 
presence of significant regions of nonequilibrium and the fact the present cases are mostly optically thin allow 
the discussion to focus on the nonequilibrium radiation (and the difference between the Boltzmann and non- 
Boltzmann models). A comparison between the measured and predicted intensity spectrum at 3. 1 cm behind the 
shock is shown in Fig. 17(a) for the 9.165 km/s case at 0.1 Torr. The measured wavelength range for this case 
is 700 - 900 nm. In addition to the EAST data and present non-Boltzmann prediction, this figure also shows the 
result of increasing the shock velocity by 1 .5% or by assuming a Boltzmann distribution of electronic states. It is 
seen that the baseline non-Boltzmann prediction under-predicts the data by nearly half. The significantly better 
prediction for the 1.5% increase in the shock velocity suggests that the 1.5% uncertainty is responsible for the 
some underprediction of the measurements. The differences seen between the non-Boltzmann and Boltzmann 
predictions confirm the nonequilibrium state of the gas at this location. 

Figure 17(b) compares the spatial variation of the integrated intensity over a sample spectral range for the 
measurement and predictions shown in Fig. 17(a). The non-Boltzmann prediction is similar in shape but lower 
in magnitude than the measurements. An important point to be made here is that measurements confirm the sup- 
pression of the significant Boltzmann-predicted emission, which is characteristic of the non-Boltzmann model 
in regions of nonequilibrium. Although the non-Boltzmann predicted nonequilibrium radiation is larger than 
the equilibrium radiation (which is approached near z = 4 cm) for these cases, the contribution of nonequilib- 
rium radiation to the wall radiative flux for CEV is small because the length of the nonequilibrium region is less 
than 2 cm along the stagnation-line, while the equilibrium region is roughly 15 cm. The Boltzmann predicted 
nonequilibrium radiation, on the other hand, is large enough in magnitude to contribute significantly to the 
wall radiative flux, even with the very small length of the nonequilibrium region. The experimental validation 
that the nonequilibrium radiation is closer to the non-Boltzmann prediction than the Boltzmann prediction is 
therefore valuable, because it confirms that the nonequilibrium radiation contribution is much less than that 
predicted by the Boltzmann model. 

4.6 Fire II, Coupled Radiation Only 

The Fire II data set from 1967 remains the touchstone for validation of entry heating simulations of radiative 
heating. [7 1-73] The vehicle had three overlay ed beryllium heatshields. A fresh surface was exposed twice 
during entry after the outermost heatshield was jettisoned as its integrity became compromised due to thermal 
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Figure 18: Sensors on Fire II vehicle, circle represents radiometer locations, x represents calorimeter locations. 




(c) Radiation intensity distribution at 1636 seconds. 




(d) Radiation intensity distribution at 1643 seconds. 


Figure 19: Total and radiative heating distributions around Fire II vehicle at trajectory points 1636 and 1643. 
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load. Consequently, radiation data was obtained without the complications of ablation products. Metallic sur- 
faces are generally considered to be strongly catalytic but the potential formation of an oxide coating introduces 
uncertainty in the proper formulation of this boundary condition. Therefore results are presented using bound- 
ing approximations for a super-catalytic condition (species mass fractions on the surface equal to free stream 
values) and for a non-catalytic condition (zero species fraction gradient at the wall). Surface temperatures were 
specified according to flight data. 

Instrument locations are presented in Fig. 18. The radiometers record radiation intensity contribution 
between 0.4 and 6.2 eV (window limits). The calorimeters record total convective and absorbed radiative 
heating. The absorbed radiative energy is calculated as a function of wavelength[74]. Simulations at two 
trajectory points are intended to illustrate capabilities at a strong thermochemical nonequilibrium condition 
(1636 seconds) and a near-equilibrium condition (1643 seconds). 

The computed heating rates for the subject points are compared to experimental data in Fig. 19. The 
total heating rates measured by the calorimeters are presented in Figs. 19 (a) and (b). The solid black line in 
these figures indicate the uncoupled total heating (radiative heating computed from converged, non-radiating 
case) assuming a super-catalytic boundary condition. The solid red line indicates the total heating with full 
coupling of radiative energy transfer. The coupling reduces the total heating by 5 to 10 % for these trajectory 
points. The dashed lines indicate the coupled, non-catalytic heating rate. The catalytic effect is most prominent 
(reduction by 1/3) at the strong nonequilibrium point at 1636 seconds. The catalytic effect is still evident at 
1643 seconds with a net reduction of 5 %. The experimental data, indicated by blue circles, is fully bounded 
by the extremes in modeling of surface catalysis. The uncoupled radiation intensity is presented with a dashed 
line and the coupled radiation intensity is presented with a solid red line in Figs. 19 (c) and (d). The uncertainty 
in the experimental data in these figures (blue circles and error bar) bounds the computed coupled results. The 
simulations indicate that radiation coupling is required to match experimental data and that the suite of physical 
models used to compute these cases are consistent with experimental data. 

4.7 Mars Return, Coupled Radiation and Ablation 

A benchmark case for modeling strongly radiating and ablating flowfields consists of a 3.05 m radius sphere 
at Voq = 15.24 km/s and - 0.000255 kg/m 3 . This case, which represents Mars return to Earth, has been 
extensively studied for over 30 years, beginning with the then state-of-the-art viscous shock layer (VSL) code 
methodology [75-77]. In Ref. [78], these previous studies were used to assess the recently implemented 
coupled radiation, coupled ablation, and free-energy minimization capability of the LAURA code [35]. The 
free-energy minimization capability was applied to provide chemical equilibrium flowfield solutions consistent 
with previous VSL solutions. The ablation rates and chemical composition assumed in the previous studies 
were applied, and the resulting stagnation line radiation, temperature, and elemental composition were shown 
to agree with the past studies. 

Since the publication of Ref. [78], three primary advancements have been made to the study of this Mars 
return example. The first of these advancements is that the assumption of chemical equilibrium has been re- 
moved, and instead the modern two-temperature thermochemical nonequilibrium model is applied. The second 
of these advancements is that the ablation rate is computed assuming steady-state ablation as part of the flow- 
field solution, and therefore the specification of the ablation rate is not required. The procedure required for 
this ablation rate computation is described by Johnston et al. [79]. Assuming the elemental mass fractions of 
injected gas are [C, H, O, N] = (0.92, 0.022, 0.049, 0.009), the predicted ablation rate at the stagnation point 
for the present Mars return case is 0.445 kg/s/m 2 . The third of these advancements is the treatment of the 
sphere’s afterbody, which allows for the treatment of the wake. For the present results, the surface of the sphere 
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Figure 20: Temperature profiles for the Mars return case with and without coupled radiation. 




Figure 21 : Radiative flux to the wall predicted with and without coupled radiation for the Mars return case. 


afterbody was assumed to be non-ablating, radiative equilibrium, and fully-catalytic to homogenous recombi- 
nation. This assumption of a non-ablating afterbody allows the influence of forebody ablation products on the 
afterbody to be studied. 

The influence of coupled radiation, without coupled ablation, on the vibrational-electronic temperature at 
various points in the flowfield is shown in Fig. 20. The radiative heating at the surface for these cases is shown 
in Fig. 21. These figures show that the radiatively cooled gas from the strongly radiating stagnation region 
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Figure 22: Mass fraction of two dominant ablation products, CO and C, for the Mars return case with coupled ablation and 
radiation. 




Figure 23: Radiative flux to the wall predicted with and without coupled ablation for the Mars return case. 


flows downstream and reduces the temperatures in the weakly radiating downstream regions of the flow. This 
non-local radiative cooling effect causes correlations such as the Tauber- Wakefield approximation[80], which 
is dependent upon the local radiative flux, to under-predict the radiative cooling effect in dowstream regions. 
This effect is apparent in Fig. 21, especially in the wake where the Tauber- Wakefield approximation predicts 
very little cooling. 

The mass fractions of CO and C predicted throughout the flowfield by assuming an ablating forebody and 
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non-abalting afterbody (as mentioned previously) are shown in Fig. 22. This figure shows that although both 
CO and C are formed as a result of the forebody ablation products, they flow into the low pressure wake and 
dominate the species composition near the body. The influence of ablation on the radiative heating is shown 
in Fig. 23. On the forebody, a reduction in the radiative heating similar to that shown previously in Ref. [78] 
is seen, while the forebody shows a surprising increase in the radiation with the introduction of ablation. This 
increase is a result of emission from the CO 4+ band system, whose contribution to the radiative flux in the 
wake is non-negligible because of the large CO concentration in that region. 


4.8 Apollo 4, Coupled Radiation and Ablation 

The Apollo 4 test case, recently reviewed by Park[8 1], provides flight data on radiative heating including effects 
of ablation products. Free stream conditions are = 10.25 km/s and = 0.000341 kg/m 3 on a 3 meter radius 
sphere. Elemental mass fractions of injected gas are [C, H, O, N] = (0.63, 0.06, 0.30, 0.01). These fractions 
were derived from data discussed by Park[81] in which quasi-steady ablation is assumed and the elemental 
fractions of char and pyrolysis gas are combined according to the fraction of their individual blowing rates. 
The wall temperature is set to 2500 K and the dimensionless blowing rate m is 0.0086 [81]. Thermochemical 
nonequilibrium model is applied. The case was initialized from a non-ablating solution in which the equilibrium 
catalytic wall boundary condition was applied. 



(a) Intensity integrated over frequency over two paths 



Figure 24: Comparisons with Apollo 4 radiometer data. 


Fig. 24 (a) plots the cumulative intensity arriving at the radiometer starting from 0.4 eV and concluding 
at 6.2 eV (the radiometer window is opaque outside this range) . The path length indicated by the solid line 
includes the shock-layer thickness immediately above the radiometer. The path length indicated by the dashed 
line includes the shock-layer thickness plus 8 cm of cavity depth (Fig. 24 (b)) which is assumed to be filled 
with gas at conditions on the surface. (There was no flow field simulation that included a cavity.) From 0.4 to 
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(a) Vibrational-electronic temperature 


(b) CO mass fraction 


Figure 25: Apollo 4 flowfield with coupled radiation and ablation. 



Afterbody 



Figure 26: Influence of coupled ablation and radiation on the forebody and afterbody convective heating for the Apollo 4 case. 


approximately 2.0 eV there is essentially no difference in the intensity reaching the radiometer over either path 
length. The intensity reaching the radiometer above 2eV becomes path length dependent, indicating that some 
of the radiation above 2eV is being absorbed by the gas in the cavity. The circular symbol with uncertainty bar 
indicates the measured cumulative intensity. It is convenient to put this symbol at the end of the window range 
at 6.2 eV but the data in fact represents the integrated intensity from 0.4 to 6.2eV. This figure shows that for 
the assumed blowing conditions the predicted cumulative intensity at the surface exceeds the measured value 
by approximately 20 %. However, if the effect of absorption over a longer path length through the cavity is 
included, then excellent agreement with measurement is obtained. 
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Figure 27: Influence of coupled ablation and radiation on the forebody radiative heating for the Apollo 4 case. 


We caution that this agreement does not constitute validation. The specified blowing rates and elemental 
mass fraction of the ablation products are estimates [81] - not computed on any first principles within this 
simulation. Still, the results confirm that for a reasonable estimate of boundary conditions there is consistency 
with measured data - providing more supporting evidence that the coupled radiation and ablation algorithms 
are correctly implemented. Further details of this simulation with more consideration of parametric variation is 
available in the literature. [78] 

Unlike the study of Ref. [78] discussed above, which approximated the Apollo stagnation region with a 3 
m hemisphere, the present results are obtained using the actual Apollo geometry at angle of attack. Elemental 
mass fractions of the char and pyrolysis gas are taken from Bartlett et al. [82] and assume mechanical silica 
removal. Steady-state ablation is assumed for the computation of the pyrolysis injection rates and wall temper- 
ature. (Recall in the previous case the surface temperature and ablation rate was specified.) This assumption 
is appropriate for the peak heating conditions considered here. The flowfield was computed assuming two- 
temperature thermochemical nonequilibrium. Similar to the Mars-return case, only the forebody is assumed to 
be ablating. The afterbody is assumed to be fully-catalytic for homogenous recombination and the temperature 
is obtained from the assumption of radiative equilibrium. 

The vibrational-electronic temperature throughout the flowfield is shown in Fig. 25 (a). This flowfield 
includes coupled radiation and ablation. The mass fractions of the dominant ablation product, CO, are shown in 
Fig. 25 (b). Note that the CO, which originates on the forebody, flows into the afterbody and result in up to 30% 
CO on the afterbody surface. The influence of this ablation product contamination is shown in Fig. 26 to result 
in a slight decrease in the convective heating on the wind- side afterbody. For lower velocity entry conditions, 
this decrease has been found to be as large as 20%. The decrease in the convective heating seen on the forebody 
in Fig. 26 is significant because this region is locally ablating. The influence of coupled radiation is shown in 
this figure to have a small impact on the convective heating. The influence of coupled ablation and radiation 
on the forebody radiative heating (the afterbody radiation is negligible) is shown in Fig 27. The influence of 
coupled radiation is shown to reduce the radiative heating by nearly 30%, while the influence of ablation reduces 
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it by only about 5%. Note that the ablation rates predicted for this case are an order-of-magnitude less than for 
the Mars return case discussed previously, in which case ablation had a significant impact on the radiation. 


4.9 Towed Ballute 

A series of wind tunnel tests were conducted at NASA Langley in Mach 6 air and CF4 to gather surface 
heating data (thermophosphor) and Schlieren video to observe unsteady motion. [16] Some of these tests were 
exploratory in nature, examining ways in which flexible, towed decelerators behaved dynamically in a ground 
based facility - accepting the fact that the Mach number and dynamic pressure are far from the flight condition. 
We discuss here a test of a hemispherically capped cylinder extending through the center of a toroid. The 
spacecraft (cylinder) diameter is 0.50 inches. The toroid diameter is 6 inches. It produces a bow shock which 
is directed just inside the inner boundary of the toroid. As the cylinder is pushed further upstream, unsteady 
motion ensues. As the cylinder is retracted, the motion appears steady (at least with regard to images captured 
in the Schlieren). The topology of the wind tunnel case is different than the flight system in that the cylindrical 
sting is not present in flight. Nevertheless, the same behavior of the recirculating flow moving forward due to 
the compression in the core of the toroid is inferred from the pulsing bow shock behind the nose of the cylinder. 

Figure 28(a) is a single frame from a Schlieren video of Run 18 in the 20-Inch Mach 6 Air Tunnel at a 
Reynolds number of 1 million per foot. The distance from the center of the hemispherical nose to the midplane 
of the toroid is 3 inches. The video frame rate is not sufficiently high to resolve a single cycle of recirculating 
flow. Instead, it captures instantaneous (to exposure time) locations of the bow shock as it responds to the 
unsteady, separated flow beneath it. The observed range of motion in this case is characterised as incipient 
unsteady response. There is a slight fluttering with small wavelength visible in the bow shock image. This 
incipient case is considered an important test because we want to confirm the ability of LAURA to identify the 
presence and extent of unsteady motion. 

A single frame from the simulation of Run 18 with LAURA is presented in Fig. 28(b) showing density 
contours. In the movie, there is significant movement of the reflected shock off the cylinder which drives the 
recirculation flow back to the nose. These structures cannot be picked up in the Schlieren. The response of the 
bow shock over the cylinder is more subdued. The fluttering of the bow shock in the CFD simulation appears 
to be somewhat stronger than observed in the experiment. The key point here is that CFD in fact picked up an 
unsteady response when one is observed in the experiment. In contrast, Run 16 (not shown here) in which the 
cylinder was retracted another inch backward, produced a flow which was steady in the Schlieren and steady in 
the CFD simulation. This agreement builds credibility for the first time that the available tools are at least able 
to discern unsteady motion when it occurs. The more complex task of validating the time accurate response is 
of less concern in that we want to avoid the unsteady motion in the flight system in the first place (or at least 
insure that the extent of unsteady motion does not compromise the flight system). 

An example of a violently unsteady response which would compromise system integrity is evident in Run 
20 in which the cylinder is pushed an additional inch forward (as compared to Run 18). In this configuration 
(Fig. 29(a)) the bow shock impinges on the toroid. A more complex wave pattern reflects off the cylinder in the 
core of the toroid, and a massively unsteady response ensues. The bow shock off the cylindrical nose moves 
from the inner to outer edge of the toroid - a situation totally unacceptable in flight from both a heating and 
system stability perspective. The unsteady simulation of Run 20 (Fig. 4.9) shows an equivalent amplitude of 
motion as the experiment. The bow shock off the spherically capped cylinder swings from the inside to outside 
edge of the trailing toroid. 
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(a) Run 18 - Schlieren (b) Simulation of Run 18 - density contours 


Figure 28: Single frame of experiment and CFD simulation for Run 18 with weakly unsteady flow (small amplitude excursions). 



(a) Run 20 - Schlieren 

Figure 29: Single frame of experiment and CFD simulation for Run 20 with strong unsteady flow (large amplitude excursions). 
The magnitude of oscillations is qualitatively captured. 
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5.0 UNCERTAINTY ANALYSES 


Numerical simulations should be expected to publish uncertainties in the same way uncertainties in experi- 
mental data are required. Simulations generally have two sources of error: (1) truncation error associated with 
spatial, temporal, and spectral discretization and (2) modeling error associated with turbulence, kinetics, trans- 
port properties, thermal relaxation, radiation, and boundary conditions. Truncation error may be quantified 
through grid convergence studies. Application of output based grid convergence (adjoint based adaptation) has 
been used to drive grid adaptation to meet user specified uncertainties in output quantities. Modeling error 
may be characterized as model form uncertainty and model parameter uncertainty. Quantification of model 
form uncertainty (e.g. algebraic turbulence model vs one- or two equation turbulence model) requires extensive 
experimental validation and the magnitudes are usually dependent on local flow conditions (mission specific, 
application specific). Quantification of uncertainty derived from model parameter uncertainties in computa- 
tional aerothermodynamic simulations has been discussed by Kleb et. al. [83, 84]. In [83] a markup format to 
input data files is described to automate uncertainty analyses in a simulation. (See Fig. 30) In [84] estimates of 
epistemic uncertainties in radiation model parameters to HARA are defined and the net uncertainty to several 
test problems are calculated — including the FIRE II simulation presented previously. The intensity between 0 
and 6 eV, above 200 nm, was measured by the Fire II radiometer and was presented in Fig. 19. The uncertainty 
analysis results in intensity of 4.8 ± 1.5 W/cm2/sr for the 1636 case, which compares well with the measured 
result of 4.9 ± 1.1 W/cm2/sr. For the 1643 case, the analysis predicts 51.1 ± 13.8 W/cm2/sr, which overlaps 
the error bounds of the measured result of 63.0 d= 14.5 W/cm2/sr. 
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Figure 30: Formats for marking up parameter uncertainties on physical property input files to LAURA. 


6.0 CONCLUDING REMARKS 

Mission requirements for entry, descent, and landing (EDL) of high mass at Mars and exploration of the outer 
planets highlight the need for new, innovative vehicle designs. Depending on the mission, simulations may 
require coupling of some combination of radiation, massive ablation, surface deflections, and jet interactions 
- sometimes in a time-accurate manner. Addressing validation of the requisite suite of physical models in 
computational aerothermodynamics (CA) remains a high priority. The need to explore innovative solutions to 
design problems has become an equally high priority. The design space for EDL systems (e.g. layout and 
thrust schedules of retro-rockets, attachment details for trailing ballutes, effect of deployable infrastructure on 
localized heating and unsteady interactions, effects of large shape change on lifting bodies) require that CA be 
brought into the design process early. These design challenges tend to introduce highly non-linear interactions 
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that are not amenable to simple engineering analysis tools. Challenges within CA discussed herein have been 
framed in the context of simulation for innovation. Simulation for innovation has three defining features: (1) 
conceptual systems can be quickly set up for simulation and optimization, (2) robust grid adaptation to a user 
specified discretization tolerance is available, and (3) credible estimates of simulated aerothermal environment 
uncertainties are published. Obstacles to any of these features impedes innovative design - be it an EDL system 
to land tens of metric tons on the Martian surface or set-up of an experiment for validation in a ground-based 
facility. 

Simulation for innovation in CA is the goal. Two algorithmic topics that must be addressed to achieve this 
goal are discussed. The first topic addresses the need to obtain accurate heating simulation on grid systems that 
are easily adapted to complex flow topologies. Tetrahedral grids provide the greatest flexibility for adaptation 
but are poorly suited for prediction of heating and shear in hypersonic flows using conventional finite- volume 
algorithms. A recent development in multi-dimensional, inviscid flux reconstruction substantially overcomes 
these limitations as demonstrated in two simple challenge problems as well as two more complex configurations 
— double cone and shuttle orbiter. The second topic addresses the infrastructure needs to support fully coupled 
ablation, radiation, and shape change in a CA simulation. A free energy minimization option enables equilib- 
rium catalytic wall and equilibrium ablation boundary conditions as well as equilibrium shock layer simulations 
in which elemental mass fractions vary due to ablation. The HARA radiation code, using the smeared-rotational 
band (SRB) method for molecular band radiation, requires only 30 to 50 seconds compute time per line of sight. 
This relative fast computation time enables coupled radiation with the flow solver for 3D configurations. 

Examples of code assessment and validation using the structured flow solver LAURA and the unstructured 
flow solver FUN3D in problems of particular interest to the CA community are presented. These cases illustrate 
important simulation capabilities but also serve to illustrate where progress is needed to achieve simulation 
for innovation. Both ground-based and flight data are computed. Simulations of shuttle orbiter heating on 
both the structured, cell-centered LAURA code formulation and on the unstructured, node-based FUN3D code 
formulation include important effects of surface catalysis. Simulations of radiative energy transfer at near- 
equilibrium and non-equilibrium conditions compare well with measured flight data (FIRE II) and ground- 
based shock tube data in EAST. Coupled massive ablation (m = 0.2) and radiation are tested on a classic Mars 
return problem with good comparisons to previous solutions using viscous shock layer methods. Comparison to 
radiometer data on Apollo 4 test flight including coupled ablation and radiation shows good agreement though 
the blowing rate in this case is prescribed. Simulation of a ground-based test of supersonic retro-propulsion 
demonstrates the difficulty and criticality of grid adaptation as judged by the plume shape and terminal shock 
shape. Simulation of unsteady flow as a function of relative position of the spacecraft and towed ballute shows 
qualitative agreement with ground-based tests. In every case, significant expert user intervention was required 
to craft a grid and/or to tune the infrastructure for multi-physics coupling. 

Finally, in an initial effort to automate the publishing of uncertainties in a given simulation, an algorithm 
is cited that allows one to compute the net uncertainty on output quantities as a function of epistemic uncer- 
tainties in physical model parameters. An example involving radiative heating is provided herein. Though this 
algorithm is unable to formally address model form uncertainties it is considered an important step in realiz- 
ing the third element of simulation for innovation — credible estimates of simulated aerothermal environment 
uncertainties are published. 

7.0 APPENDIX A - MULTI DIMENSIONAL RECONSTRUCTION 

The multi-dimensional reconstruction algorithm used to enable a tetrahedral grid system for surface heating is 
provided here for the convenience of the reader. More details are available in the original documentation. [31] 


RTO-EN-AVT-1 86 


13- 35 



Challenges to Computational Aerothermodynamic Simulation 
and Validation for Planetary Entry Vehicle Analysis 



ORGANIZATION 


The new algorithm is implemented within the flow solver FUN3D - a node based, fully unstructured, finite- 
volume solver of the Euler and Navier-Stokes equations. [85] The method uses a Green-Gauss formulation of 
gradients at the nodes for multi-dimensional, inviscid flux reconstruction. Viscous gradients are also computed 
from a Green-Gauss formulation. A suite of modules includes all of the gas physics models in LAURA and 
VULCAN[86] for thermodynamics, transport properties, chemical kinetics, and thermal relaxation. The base- 
line inviscid flux reconstruction algorithms utilize quasi-one-dimensional (edge-based) reconstruction using 
Roe’s averaging with Yee’s symmetric total variation diminishing algorithm adapted for unstructured grids. 

A triangular grid system is presented in Fig. 31. The two-dimensional reconstruction algorithm is developed 
relative to this figure. This derivation is repeated from the original source [29] for the convenience of the reader 
with special emphasis placed on updates to the algorithm. 

Assume a node based scheme as illustrated in Fig. 31 with (x, y) coordinates of all nodes (1,2,3, . . . ) given 
and dependent variables q at these same nodes to be determined by solution of conservation laws. The conser- 
vation laws are discretized relative to the dual control volume, represented by dashed lines in Fig. 31. The dual 
volume around node 3 passes through the centroids of all elements surrounding the node (A, L>, C, 19, E , F ) 
and the midpoints of all edges separating these elements. 



Figure 31 : Schematic of node-based, unstructured grid system showing element centroids (uppercase letters), edge midpoints 
(lowercase letters), nodes (numerals), the dual volume (dashed lines), and the principal directions (dotted lines with arrowhead). 


7.0.1 Reconstruction Overview - Principal (Rotated) Direction 

Consider a conventional, quasi-one-dimensional, first-order reconstruction of inviscid flux across face AcB 
separating the dual volumes around nodes 1 and 3. This flux is defined as a function of information at nodes 1 
and 3 of edge c: 

f AcB = 1 [f3 + fi - R-c|A c |R~ 1 (q 3 - qi)] (20) 

Here R c and A c are the matrix of eigenvectors and diagonal matrix of eigenvalues for the Jacobian of f with 
respect to q. Inviscid fluxes are then defined using a loop over edges and the reconstruction direction is fully 
determined by the grid. 
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The basic idea in multi-dimensional reconstruction is that: (1) reconstruction directions should be based on 
local flow characteristics and not be constrained by the grid; and (2) the flux components for any face will utilize 
a stencil employing all nodes of a surrounding element. Conventional, node-based schemes loop over edges 
when using quasi-one-dimensional reconstruction. In contrast, the new multi-dimensional reconstruction loops 
over elements. Consider the two-dimensional element defined by nodes (1, 2, 3) with centroid A in Fig. 31. 
The computation of inviscid flux in this element employs the same infrastructure as the viscous flux. 

A principal direction in the element, x', is defined by the direction of Vp, computed using the same Green 
Gauss evaluation of the derivatives that are already used to compute the viscous terms. If density is locally 
constant then the principal direction aligns with the local velocity computed as the average of nodal values. The 
inviscid flux within the element will be computed along and orthogonal to the principal flow direction. 

The flux components are defined using nodal weights based on the Green-Gauss formulation for derivatives 
in the component directions. Thus, 


df a 

dx f 


1 

Qa 


^ ^ ^k^k^k,x' 

k= faces 


( 21 ) 


where is the average value of f across all nodes defining surface k of element A. With reference to Fig. 31 
this derivative can be written 


dl, 


x' ,A 
dx f 


— a 3,x'{h,x' + + a l : x'(^3,x' + ^2,x') + a 2,x'(^l,x' + ?3,x') 


where 


°tk,x' 


^k^k^x 

2VL A 


and k is a face index identified by the number of the opposite node. The sum J2k=faces 
closed element. It is convenient to reorder and scale the coefficients as follows. 

df x ', A 


dx' 


— ( a 2,x' + a 3,x')fl + ( a 3,x' + a l,x')^2 + (Oi i^ x ' + ^2,x')^3 


+ ^2,x'h + 03,x f k] 


Ax' 


Note that Ax' can be interpreted as a distance across the element in x' direction. It is defined 

2 


Ax' 


— ^ ^ (|^n— l,x' T" | 


n=nodes 


(22) 

(23) 
0 for any 

(24) 

(25) 

(26) 


where a cyclic indexing is assumed. Note that the original formulation [29] is corrected here with a factor 
2 in the numerator which provides the distance between virtual nodes. The reordering and scaling yield the 
following relations for /3 n ^ x r. 

@n,x' ~ n—l,x ' &n- 1-1,#') (27) 

= 0 (28) 

n 

£|M = 2 (29) 

n 

Two options for computing i x > are developed. Both mirror the baseline quasi-one-dimension reconstruction 
using Roe’s averaging across nodes. The first option averages dependent variables at two virtual nodes used in 
the reconstruction. The second option may be considered a weighted average of a rotated difference scheme on 
existing edges. No virtual nodes are created. 
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7.0.2 Option 1 - Virtual Node Averaging 


A right and left virtual state for computing i x < are computed as follows. 


Q.R,x' ^ ^ (l/^n^'l “b fin,xfityji 

(30) 

n=nodes 


O.L,x' = ^ ^ (lAi,,#'! — fin,x')Qn 

(31) 

n=nodes 


fx' = \ [f(q)L,a!' + ~ dq x ',Um) 

(32) 


where the left and right states of f are functions of the left and right states of q, respectively. The matrices 
R x / and A um, x ' are computed as functions of the Roe’s average of the right and left states. (Forming q R iX r 
and q^/ from averages of l/pi, u, v, w, 1/T have provided the most robust simulations as judged by level of 
temperature undershoots across strong bow shocks.) 


dc{x' 


F ^x' ( qr,x' 

FF ^ 7 ^ ^ fin,x'Qn 
n=nodes 


— Ax Ft x / ^ ^ (|/ 7 n,a?'| “b fin,x' Qn,GG 

n=nodes 

dc\L,x' Ax R x / ^ ^ (IAt,,#'! fin,xfi^ 7 Qn,GG 

n=nodes 


dcix'jim = $ minmod 


‘^dc^R xG ^dc^x' i ^dc^R xG 2 (dc^R^x' ~b dciL xfi 


where <f> is an auxiliary limiter (0 < <J> < 1) defined in a later section. 


(33) 


(34) 

(35) 

(36) 


7.0.3 Option 2 - Weighted Average of Edges to Principal Node 

Select the largest of \fi n ,x' | to identify the principal node for construction of f x ' and express it as a function of 
the remaining coefficients. Assume node 3 is the principal node. Then 

df x ' = fil,x'fl + fi2,x'f2 + fi^x'h 

— fil,x'fl + fi2,x'f2 — (fil,x f + fi2,x')f3 (37) 

= fii,x' (fi ~ h) + fi2,x f (/2 - h) 

The reconstructed flux in the direction x', f^/, is computed as a weighted average of surrounding edges. Fur- 
thermore, it is noted that if any of the surrounding edges are parallel to x' then the weight of that edge will 
equal 1 and the weight of the other edges will equal zero. 

f x ' = \fil,x'\h-3,x' + \ fi2,x'\h-3,x' (38) 


fl-3,x' 

h-3,x' 


1 

2 

1 

2 


^l,x' “b ^ 3 ,x' s iS n (/ 5 l,a; / )Fti_ 3 ^/ 1 ^lim,l— 3 ,x' | (^ 0 . 1 — 3 ,^' 

F 2 ,x ' “b h,x' — si&^-(fi 2 ,x')^ j 2 — 3 ,x' I A/zm, 2 — 3 ,x' \ (d(\ 2 — 3 ,x' 


dQl—3,x',lim) 

dc[2—3,x',lim) 


(39) 

(40) 
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d£[i— 3^' 

— Ri-3,x'(qi • 

- qs) 


(41) 

rfOL 2-3, x' 

= R 2 -3,x.'(q 2 ■ 

-q3) 


(42) 

dQl —3,x',lim 

= $ minmod 

2,dc[i x ' * 2dq3 x i. 

, ^(dqi,*' + <*03,*') 

(43) 

dc\2 —3,x',lim 

= $ minmod 

2dq2,a; , > ^dci2-3,x'i 2dq3 ja; '. 

, 1 (dq 2 ,x' + <*03,*') 

(44) 

dc{n,x' 

— Ax 


(45) 


Note again that is the same auxiliary limiter used in Option 1. 
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